# Author: Patrick Cunha Silva
# Project: Janusz, Cunha, and Junqueira
# Goal: SI for Afro-Brazilian Analysis
# First Version: Jun 21, 2023
# This Version: Jun 4, 2024

rm(list = ls())

# Read genFunctions
source(here("R-Files", "genFunctions.R"))

# Load data
load(file = here("Data", "blackRDD_final.RData")) # main data (RD Models)
load(file = here("Data", "candRace.RData")) # data for racial composition of the 2016 Election
load(here("Data","df_spatial_nb_std_race_final.RData")) # for spatial analysis


# Session ran on 06/04/2024
# > sessionInfo()
# R version 4.3.0 (2023-04-21)
# Platform: aarch64-apple-darwin22.4.0 (64-bit)
# Running under: macOS 14.5
# 
# Matrix products: default
# BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
# LAPACK: /opt/homebrew/Cellar/r/4.3.0_1/lib/R/lib/libRlapack.dylib;  LAPACK version 3.11.0
# 
# locale:
#   [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# time zone: America/Chicago
# tzcode source: internal
# 
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] stargazer_5.2.3 rdpower_2.2     rddensity_2.4   rdrobust_2.1.1  here_1.0.1     
# 
# loaded via a namespace (and not attached):
#   [1] vctrs_0.6.5      cli_3.6.1        rlang_1.1.1      generics_0.1.3   lpdensity_2.4    glue_1.6.2       colorspace_2.1-0 rprojroot_2.0.3 
# [9] scales_1.2.1     fansi_1.0.4      grid_4.3.0       munsell_0.5.0    tibble_3.2.1     MASS_7.3-58.4    lifecycle_1.0.3  compiler_4.3.0  
# [17] dplyr_1.1.2      pkgconfig_2.0.3  rstudioapi_0.14  R6_2.5.1         tidyselect_1.2.0 utf8_1.2.3       pillar_1.9.0     magrittr_2.0.3  
# [25] tools_4.3.0      gtable_0.3.3     ggplot2_3.4.2   



#######################################################
########### Part 1 - Descriptive Statistics ###########
#######################################################

# Racial Composition of the 2016 Mayoral Election (winner and runner-up)
table(winner = candRace$DS_COR_RACA_WINNER, runnerup = candRace$DS_COR_RACA_RUNNERUP)

# variables to summary stats
vars <- c('run_var_black', 'totalBlacks', 'totalBlacks_ratio', 'totalBlacks_pct',
          'totalBlacks_noexp', 'totalBlacks_noexp_ratio', 'totalBlacks_noexp_pct',
          'totalBlacks_noexp_frommp_pct', 'totalBlacks_noexp_frommp_ratio', 'totalBlacks_noexp_frommp_pct')

stargazer(blackRDD[, vars],
          covariate.labels = c('Running Variable', 
                               '# of Afro-Brazilian Candidates',
                               'Ratio of Afro-Brazilian Candidates Per Seat',
                               '% of Afro-Brazilian Candidates Per Seat',
                               '# of Afro-Brazilian Candidates (No Political Experience)',
                               'Ratio of Afro-Brazilian Candidates Per Seat (No Political Experience)',
                               '% of Afro-Brazilian Candidates Per Seat (No Political Experience)',                               
                               '# of Afro-Brazilian Candidates (No Political Experience from Mayoral Party)',
                               'Ratio of Afro-Brazilian Candidates Per Seat (No Political Experience from Mayoral Party)',
                               '% of Afro-Brazilian Candidates Per Seat (No Political Experience from Mayoral Party)'                               
                               ),
          title = 'Descriptive Statistics for the Afro-Brazilian Candidate Emergence Analysis',
          type = 'text'
          )

# see file 7_SI_DescriptivePlot.R for figures

##############################################################
########### Part 2 - Models With Control Variables ###########
##############################################################

# control variables
control.vars <- c('woman', 'times.as.candidate', 'age', 'times.as.winner',
                  'incumbent', 'married', 'college', 'merchant',
                  'politician', 'servant', 'business', 'doctor',
                  'administrator', 'farmer',
                  'black_pop_pct', 'urban_pct', 'ln_gdp18_pc', 'female_pct', 'population2018')

############### No Restrictions

# Black Candidates (Ratio)
m1a <- rdrobust(y = blackRDD$totalBlacks,
                x = blackRDD$run_var_black,
                covs = blackRDD[, control.vars],
                all = TRUE)
summary(m1a) # Positive, Null


# Black Candidates (Ratio)
m1b <- rdrobust(y = blackRDD$totalBlacks_ratio,
                x = blackRDD$run_var_black,
                covs = blackRDD[, control.vars],                
                all = TRUE)
summary(m1b) # Positive, Null


# Black Candidates (pct)
m1c <- rdrobust(y = blackRDD$totalBlacks_pct,
                x = blackRDD$run_var_black,
                covs = blackRDD[, control.vars],                  
                all = TRUE)
summary(m1c) # Negative, Null


############### No Experience

# Black Candidates
m2a  <- rdrobust(y = blackRDD$totalBlacks_noexp,
                 x = blackRDD$run_var_black,
                 covs = blackRDD[, control.vars],                  
                 all = TRUE)
summary(m2a) # Negative, Null

# Black Candidates (ratio)
m2b  <- rdrobust(y = blackRDD$totalBlacks_noexp_ratio,
                 x = blackRDD$run_var_black,
                 covs = blackRDD[, control.vars],                  
                 all = TRUE)
summary(m2b) # Negative, Null

# Black Candidates (pct)
m2c  <- rdrobust(y = blackRDD$totalBlacks_noexp_pct,
                 x = blackRDD$run_var_black,
                 covs = blackRDD[, control.vars],                  
                 all = TRUE)
summary(m2c) # Positive, Null

############### No Experience From Mayoral Party

# Black Candidates
m3a  <- rdrobust(y = blackRDD$totalBlacks_noexp_frommp,
                 x = blackRDD$run_var_black,
                 covs = blackRDD[, control.vars],                  
                 all = TRUE)
summary(m3a) # Negative, p<0.05

# Black Candidates (ratio)
m3b  <- rdrobust(y = blackRDD$totalBlacks_noexp_frommp_ratio,
                 x = blackRDD$run_var_black,
                 covs = blackRDD[, control.vars],                  
                 all = TRUE)
summary(m3b) # Negative, p<0.05

# Black Candidates (pct)
m3c  <- rdrobust(y = blackRDD$totalBlacks_noexp_frommp_pct,
                 x = blackRDD$run_var_black,
                 covs = blackRDD[, control.vars],                  
                 all = TRUE)
summary(m3c) # Negative, Null

#####################################################
########### Part 3 - Sensitivity Analysis ###########
#####################################################

############### No Restrictions
pdf(here('plots_tables', 'Sens.BlackCandidates.pdf'), height = 8, width = 8)
sensitivity(mydata = blackRDD, 
            DV = 'totalBlacks', 
            run_var = 'run_var_black',
            ylabel = '# of Afro-Brazilian Candidates',
            ylb = -60, yub = 60)
sensitivity(mydata = blackRDD, 
            DV = 'totalBlacks_ratio', 
            run_var = 'run_var_black',
            ylabel = 'Ratio of Afro-Brazilian Candidates Per Seat',
            ylb = -10, yub = 10)
sensitivity(mydata = blackRDD, 
            DV = 'totalBlacks_pct', 
            run_var = 'run_var_black',
            ylabel = '% of Afro-Brazilian Candidates',
            ylb = -60, yub = 60)
dev.off()


############### No Experience
pdf(here('plots_tables', 'Sens.BlackCandidates_NoExp.pdf'), height = 8, width = 8)
sensitivity(mydata = blackRDD, 
            DV = 'totalBlacks_noexp', 
            run_var = 'run_var_black',
            ylabel = '# of Afro-Brazilian Candidates',
            ylb = -60, yub = 60)
sensitivity(mydata = blackRDD, 
            DV = 'totalBlacks_noexp_ratio', 
            run_var = 'run_var_black',
            ylabel = 'Ratio of Afro-Brazilian Candidates Per Seat',
            ylb = -10, yub = 10)
sensitivity(mydata = blackRDD, 
            DV = 'totalBlacks_noexp_pct', 
            run_var = 'run_var_black',
            ylabel = '% of Afro-Brazilian Candidates',
            ylb = -60, yub = 60)
dev.off()

############### No Experience From Mayoral Party
pdf(here('plots_tables', 'Sens.BlackCandidates_NoExpMayoral.pdf'), height = 8, width = 8)
sensitivity(mydata = blackRDD, 
            DV = 'totalBlacks_noexp_frommp', 
            run_var = 'run_var_black',
            ylabel = '# of Afro-Brazilian Candidates',
            ylb = -15, yub = 15)
sensitivity(mydata = blackRDD, 
            DV = 'totalBlacks_noexp_frommp_ratio', 
            run_var = 'run_var_black',
            ylabel = 'Ratio of Afro-Brazilian Candidates Per Seat',
            ylb = -1, yub = 1)
sensitivity(mydata = blackRDD, 
            DV = 'totalBlacks_noexp_frommp_pct', 
            run_var = 'run_var_black',
            ylabel = '% of Afro-Brazilian Candidates',
            ylb = -60, yub = 60)
dev.off()


#############################################################################
########### Part 4 - Models fo Cities with >= 50% Afro-Brazilians ###########
#############################################################################


############### No Restrictions

# Black Candidates (Ratio)
m1a <- rdrobust(y = blackRDD[blackRDD$black_pop_pct >= 50, ]$totalBlacks,
                x = blackRDD[blackRDD$black_pop_pct >= 50, ]$run_var_black,
                all = TRUE)
summary(m1a) # Positive, Null


# Black Candidates (Ratio)
m1b <- rdrobust(y = blackRDD[blackRDD$black_pop_pct >= 50, ]$totalBlacks_ratio,
                x = blackRDD[blackRDD$black_pop_pct >= 50, ]$run_var_black,
                all = TRUE)
summary(m1b) # Positive, Null


# Black Candidates (pct)
m1c <- rdrobust(y = blackRDD[blackRDD$black_pop_pct >= 50, ]$totalBlacks_pct,
                x = blackRDD[blackRDD$black_pop_pct >= 50, ]$run_var_black,
                all = TRUE)
summary(m1c) # Negative, Null


############### No Experience

# Black Candidates
m2a  <- rdrobust(y = blackRDD[blackRDD$black_pop_pct >= 50, ]$totalBlacks_noexp,
                 x = blackRDD[blackRDD$black_pop_pct >= 50, ]$run_var_black,
                 all = TRUE)
summary(m2a) # Positive, Null

# Black Candidates (ratio)
m2b  <- rdrobust(y = blackRDD[blackRDD$black_pop_pct >= 50, ]$totalBlacks_noexp_ratio,
                 x = blackRDD[blackRDD$black_pop_pct >= 50, ]$run_var_black,
                 all = TRUE)
summary(m2b) # Positive, Null

# Black Candidates (pct)
m2c  <- rdrobust(y = blackRDD[blackRDD$black_pop_pct >= 50, ]$totalBlacks_noexp_pct,
                 x = blackRDD[blackRDD$black_pop_pct >= 50, ]$run_var_black,
                 all = TRUE)
summary(m2c) # Positive, Null

############### No Experience From Mayoral Party

# Black Candidates
m3a  <- rdrobust(y = blackRDD[blackRDD$black_pop_pct >= 50, ]$totalBlacks_noexp_frommp,
                 x = blackRDD[blackRDD$black_pop_pct >= 50, ]$run_var_black,
                 all = TRUE)
summary(m3a) # Negative, Null

# Black Candidates (ratio)
m3b  <- rdrobust(y = blackRDD[blackRDD$black_pop_pct >= 50, ]$totalBlacks_noexp_frommp_ratio,
                 x = blackRDD[blackRDD$black_pop_pct >= 50, ]$run_var_black,
                 all = TRUE)
summary(m3b) # Negative, p<0.05

# Black Candidates (pct)
m3c  <- rdrobust(y = blackRDD[blackRDD$black_pop_pct >= 50, ]$totalBlacks_noexp_frommp_pct,
                 x = blackRDD[blackRDD$black_pop_pct >= 50, ]$run_var_black,
                 all = TRUE)
summary(m3c) # Negative, Null



#################################################################
########### Part 5 - Models For Different Polynomials ###########
#################################################################

# Change the Polynomial value
polynomial <- 2

############### No Restrictions

# Black Candidates (Ratio)
m1a <- rdrobust(y = blackRDD$totalBlacks,
                x = blackRDD$run_var_black,
                p = polynomial,
                all = TRUE)
summary(m1a) 


# Black Candidates (Ratio)
m1b <- rdrobust(y = blackRDD$totalBlacks_ratio,
                x = blackRDD$run_var_black,
                p = polynomial,                
                all = TRUE)
summary(m1b) 


# Black Candidates (pct)
m1c <- rdrobust(y = blackRDD$totalBlacks_pct,
                x = blackRDD$run_var_black,
                p = polynomial,                
                all = TRUE)
summary(m1c) 


############### No Experience

# Black Candidates
m2a  <- rdrobust(y = blackRDD$totalBlacks_noexp,
                 x = blackRDD$run_var_black,
                 p = polynomial,                 
                 all = TRUE)
summary(m2a) 

# Black Candidates (ratio)
m2b  <- rdrobust(y = blackRDD$totalBlacks_noexp_ratio,
                 x = blackRDD$run_var_black,
                 p = polynomial,                 
                 all = TRUE)
summary(m2b) 

# Black Candidates (pct)
m2c  <- rdrobust(y = blackRDD$totalBlacks_noexp_pct,
                 x = blackRDD$run_var_black,
                 p = polynomial,                 
                 all = TRUE)
summary(m2c)

############### No Experience From Mayoral Party

# Black Candidates
m3a  <- rdrobust(y = blackRDD$totalBlacks_noexp_frommp,
                 x = blackRDD$run_var_black,
                 p = polynomial,               
                 all = TRUE)
summary(m3a) 

# Black Candidates (ratio)
m3b  <- rdrobust(y = blackRDD$totalBlacks_noexp_frommp_ratio,
                 x = blackRDD$run_var_black,
                 p = polynomial,                 
                 all = TRUE)
summary(m3b) 

# Black Candidates (pct)
m3c  <- rdrobust(y = blackRDD$totalBlacks_noexp_frommp_pct,
                 x = blackRDD$run_var_black,
                 p = polynomial,                
                 all = TRUE)
summary(m3c) 

###############################################
########### Part 6 - Power Analysis ###########
###############################################

############### No Restrictions
p1a <- rdsampsi(data = na.omit(blackRDD[, c('totalBlacks', 'run_var_black')]), 
         tau = sd(blackRDD$totalBlacks, na.rm = TRUE)/2, cutoff = 0)
p1b <- rdsampsi(data = na.omit(blackRDD[, c('totalBlacks_ratio', 'run_var_black')]), 
               tau = sd(blackRDD$totalBlacks_ratio, na.rm = TRUE)/2, cutoff = 0)
p1c <- rdsampsi(data = na.omit(blackRDD[, c('totalBlacks_pct', 'run_var_black')]), 
               tau = sd(blackRDD$totalBlacks_pct, na.rm = TRUE)/2, cutoff = 0)

############### No Restrictions
p2a <- rdsampsi(data = na.omit(blackRDD[, c('totalBlacks_noexp', 'run_var_black')]), 
                tau = sd(blackRDD$totalBlacks_noexp, na.rm = TRUE)/2, cutoff = 0)
p2b <- rdsampsi(data = na.omit(blackRDD[, c('totalBlacks_noexp_ratio', 'run_var_black')]), 
                tau = sd(blackRDD$totalBlacks_noexp_ratio, na.rm = TRUE)/2, cutoff = 0)
p2c <- rdsampsi(data = na.omit(blackRDD[, c('totalBlacks_noexp_pct', 'run_var_black')]), 
                tau = sd(blackRDD$totalBlacks_noexp_pct, na.rm = TRUE)/2, cutoff = 0)

############### No Restrictions
p3a <- rdsampsi(data = na.omit(blackRDD[, c('totalBlacks_noexp_frommp', 'run_var_black')]), 
                tau = sd(blackRDD$totalBlacks_noexp_frommp, na.rm = TRUE)/2, cutoff = 0)
p3b <- rdsampsi(data = na.omit(blackRDD[, c('totalBlacks_noexp_frommp_ratio', 'run_var_black')]), 
                tau = sd(blackRDD$totalBlacks_noexp_frommp_ratio, na.rm = TRUE)/2, cutoff = 0)
p3c <- rdsampsi(data = na.omit(blackRDD[, c('totalBlacks_noexp_frommp_pct', 'run_var_black')]), 
                tau = sd(blackRDD$totalBlacks_noexp_frommp_pct, na.rm = TRUE)/2, cutoff = 0)

#################################################
########### Part 7 - Assumption Tests ###########
#################################################

#-------- Density Test ------#
# Directory to save plots
summary(t1 <- rddensity(blackRDD$run_var_black, c = 0)) # No manipulation
# Result: No sorting
# Density Plot 
pdf(here('plots_tables', 'Blacks_density.pdf'), height = 8, width = 12)
par(mar = c(6, 6, 2, 2))
hist(blackRDD$run_var_black, breaks = 75, # at the bw
     main = "", xlab = "Difference in Vote Share", 
     xlim = c(-50, 50),
     ylim = c(0, 80),
     cex.lab = 1.75, 
     cex.axis = 1.5,
     col = 'grey90')
abline(v = 0, col = "black", lwd = 8)
text(x = 30, y = 75, labels = paste("p-value = ", round(t1$test$p_jk, 3)), cex = 1.75)
dev.off()

# Balance Tests
summary(rdrobust(y = blackRDD$woman,
        x = blackRDD$run_var_black,
        all = TRUE))
summary(rdrobust(y = blackRDD$college,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$married,
        x = blackRDD$run_var_black,
        all = TRUE))
summary(rdrobust(y = blackRDD$age,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$merchant,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$servant,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$administrator,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$farmer,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$business,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$doctor,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$politician,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$incumbent,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$times.as.candidate,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$times.as.winner,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$black_pop_pct,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$female_pct,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$urban_pct,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$ln_gdp18_pc,
               x = blackRDD$run_var_black,
               all = TRUE))
summary(rdrobust(y = blackRDD$population2018,
                 x = blackRDD$run_var_black,
                 all = TRUE))


#############################################################
########### Part 7 - Models For Pretos and Pardos ###########
#############################################################


############### Pretos
# Pretos Candidates (Ratio)
m1a <- rdrobust(y = blackRDD$totalPreto[blackRDD$DS_COR_RACA_RUNNERUP == 'PRETA' | blackRDD$DS_COR_RACA_WINNER == 'PRETA'],
                x = blackRDD$run_var_black[blackRDD$DS_COR_RACA_RUNNERUP == 'PRETA' | blackRDD$DS_COR_RACA_WINNER == 'PRETA'],
                all = TRUE)
summary(m1a) # Negative,Null


# Pretos Candidates (Ratio)
m1b <- rdrobust(y = blackRDD$totalPreto_ratio[blackRDD$DS_COR_RACA_RUNNERUP == 'PRETA' | blackRDD$DS_COR_RACA_WINNER == 'PRETA'],
                x = blackRDD$run_var_black[blackRDD$DS_COR_RACA_RUNNERUP == 'PRETA' | blackRDD$DS_COR_RACA_WINNER == 'PRETA'],
                all = TRUE)
summary(m1b) # Negative, Null


# Pretos Candidates (pct)
m1c <- rdrobust(y = blackRDD$totalPreto_pct[blackRDD$DS_COR_RACA_RUNNERUP == 'PRETA' | blackRDD$DS_COR_RACA_WINNER == 'PRETA'],
                x = blackRDD$run_var_black[blackRDD$DS_COR_RACA_RUNNERUP == 'PRETA' | blackRDD$DS_COR_RACA_WINNER == 'PRETA'],
                all = TRUE)
summary(m1c) # Negative, p<0.05

############### Browns
# Browns Candidates (Ratio)
m1a <- rdrobust(y = blackRDD$totalPardo[blackRDD$DS_COR_RACA_RUNNERUP == 'PARDA' | blackRDD$DS_COR_RACA_WINNER == 'PARDA'],
                x = blackRDD$run_var_black[blackRDD$DS_COR_RACA_RUNNERUP == 'PARDA' | blackRDD$DS_COR_RACA_WINNER == 'PARDA'],
                all = TRUE)
summary(m1a) # Positive, Null


# Browns Candidates (Ratio)
m1b <- rdrobust(y = blackRDD$totalPardo_ratio[blackRDD$DS_COR_RACA_RUNNERUP == 'PARDA' | blackRDD$DS_COR_RACA_WINNER == 'PARDA'],
                x = blackRDD$run_var_black[blackRDD$DS_COR_RACA_RUNNERUP == 'PARDA' | blackRDD$DS_COR_RACA_WINNER == 'PARDA'],
                all = TRUE)
summary(m1b) # Positive, Null


# Browns Candidates (pct)
m1c <- rdrobust(y = blackRDD$totalPardo_pct[blackRDD$DS_COR_RACA_RUNNERUP == 'PARDA' | blackRDD$DS_COR_RACA_WINNER == 'PARDA'],
                x = blackRDD$run_var_black[blackRDD$DS_COR_RACA_RUNNERUP == 'PARDA' | blackRDD$DS_COR_RACA_WINNER == 'PARDA'],
                all = TRUE)
summary(m1c) # Negative, Null

################################################################################
### Part 10 - Candidates with Experience V. Candidates Without Experience ######
################################################################################

#### No Experience
# Black Candidates
mNo1  <- rdrobust(y = blackRDD$totalBlacks_noexp,
                 x = blackRDD$run_var_black,
                 all = TRUE)
summary(mNo1) # Positive, Null

# Black Candidates (ratio)
mNo2  <- rdrobust(y = blackRDD$totalBlacks_noexp_ratio,
                 x = blackRDD$run_var_black,
                 all = TRUE)
summary(mNo2) # Positive, Null

# Black Candidates (pct)
mNo3  <- rdrobust(y = blackRDD$totalBlacks_noexp_pct,
                 x = blackRDD$run_var_black,
                 all = TRUE)
summary(mNo3) # Negative, Null

#### With Experience
# Black Candidates
mW1  <- rdrobust(y = blackRDD$totalBlacks_exp,
                  x = blackRDD$run_var_black,
                  all = TRUE)
summary(mW1) # Positive, Null

# Black Candidates (ratio)
mW2  <- rdrobust(y = blackRDD$totalBlacks_exp_ratio,
                  x = blackRDD$run_var_black,
                  all = TRUE)
summary(mW2) # Positive, Null

# Black Candidates (pct)
mW3  <- rdrobust(y = blackRDD$totalBlacks_exp_pct,
                  x = blackRDD$run_var_black,
                  all = TRUE)
summary(mW3) # Positive, Null


# Difference in Coefficients
zDif(beta1 = mNo1$Estimate[1, 2], se1 = mNo1$Estimate[1, 4], beta2 = mW1$Estimate[1, 2], se2 = mW1$Estimate[1, 4])
zDif(beta1 = mNo2$Estimate[1, 2], se1 = mNo2$Estimate[1, 4], beta2 = mW2$Estimate[1, 2], se2 = mW2$Estimate[1, 4])
zDif(beta1 = mNo3$Estimate[1, 2], se1 = mNo3$Estimate[1, 4], beta2 = mW3$Estimate[1, 2], se2 = mW3$Estimate[1, 4])

##################################
### Part 9 - Spatial Models ######
##################################

# Immediate neighbors #
# All black candidates #

spec = formula(n_black_cands_council ~ spatial_black_mayor +
                 lag_black_elected_mayor+
                 lag_n_black_elected_council + 
                 lag_n_black_cands_council +
                 lag_perc_votes_black_council +
                 n_vagas_council +
                 n_cands_council +
                 lag_perc_votes_left_council +
                 perc_literate +
                 perc_female +
                 perc_urban +
                 perc_white +
                 log_pop +
                 pib_pc +
                 as.factor(uf))

m_poisson_black_all = glm(spec, family="poisson", data=df_spatial_nb_std_race)
m_poisson_black_all_cl = coeftest(m_poisson_black_all, vcov = vcovCL, cluster = ~uf) # standard errors clustered by state

m_neg_bin_black_all = glm.nb(spec, data=df_spatial_nb_std_race)
m_neg_bin_black_all_cl = coeftest(m_neg_bin_black_all, vcov = vcovCL, cluster = ~uf) # standard errors clustered by state

spec_ols = formula(log_n_black_cands_council ~ spatial_black_mayor +
                     lag_black_elected_mayor+
                     lag_n_black_elected_council + 
                     lag_n_black_cands_council +
                     lag_perc_votes_black_council +
                     n_vagas_council +
                     n_cands_council +
                     lag_perc_votes_left_council +
                     perc_literate +
                     perc_female +
                     perc_urban +
                     perc_white +
                     log_pop +
                     pib_pc +
                     as.factor(uf))

m_ols_black_all = lm(spec_ols, data=df_spatial_nb_std_race)
m_ols_black_all_cl = coeftest(m_ols_black_all, vcov = vcovCL, cluster = ~uf) # standard errors clustered by state

texreg(list(m_poisson_black_all_cl,m_neg_bin_black_all_cl,m_ols_black_all_cl))

# Immediate neighbors #
# All black with no experience #

spec = formula(n_cands_black_no_exp_council ~ spatial_black_mayor +
                 lag_black_elected_mayor+
                 lag_n_black_elected_council + 
                 lag_n_black_cands_council +
                 lag_perc_votes_black_council +
                 n_vagas_council +
                 n_cands_council +
                 lag_perc_votes_left_council +
                 perc_literate +
                 perc_female +
                 perc_urban +
                 perc_white +
                 log_pop +
                 pib_pc +
                 as.factor(uf))

m_poisson_black_noexp = glm(spec, family="poisson", data=df_spatial_nb_std_race)
m_poisson_black_noexp_cl = coeftest(m_poisson_black_noexp, vcov = vcovCL, cluster = ~uf) # standard errors clustered by state

m_neg_bin_black_noexp = glm.nb(spec, data=df_spatial_nb_std_race)
m_neg_bin_black_noexp_cl = coeftest(m_neg_bin_black_noexp, vcov = vcovCL, cluster = ~uf) # standard errors clustered by state

spec_ols = formula(log_n_cands_black_no_exp_council ~ spatial_black_mayor +
                     lag_black_elected_mayor+
                     lag_n_black_elected_council + 
                     lag_n_black_cands_council +
                     lag_perc_votes_black_council +
                     n_vagas_council +
                     n_cands_council +
                     lag_perc_votes_left_council +
                     perc_literate +
                     perc_female +
                     perc_urban +
                     perc_white +
                     log_pop +
                     pib_pc +
                     as.factor(uf))

m_ols_black_noexp = lm(spec_ols, data=df_spatial_nb_std_race)
m_ols_black_noexp_cl = coeftest(m_ols_black_noexp, vcov = vcovCL, cluster = ~uf)

texreg(list(m_poisson_black_noexp_cl,m_neg_bin_black_noexp_cl,m_ols_black_noexp_cl))

